Measurement fault detection

ABSTRACT

The present invention discloses a method and an apparatus for improving measurement fault detection in a sequential measurement processing estimator, and is particularly applied to Global Positioning Receivers.

BACKGROUND OF THE INVENTION

[0001] 1. Field of the Invention

[0002] The present invention relates in general to the use ofvariable-gain filters, such as Kalman filters to estimate the state of asystem, and in particular to improved measurement fault detection insuch filters.

[0003] 2. Description of the Related Art

[0004] The discrete time Kalman filter was introduced by Rudolf Kalmanin a “A New Approach to Linear Filtering and Prediction Problems,” inthe Transactions of the ASME, Series D: Journal of Basic Engineering,Vol. 82, pp. 3545 in 1960. The Kalman filter (KF) produces the stateestimate with the smallest mean-square error for a linear time-invariantsystem with white, Gaussian state and measurement noise. The Kalmanfilter extended the concept of least-squares estimation to dynamicsystems. The concept of least-squares estimation dates back to Gauss, atleast. The Kalman filter has found innumerable applications innavigation, tracking, process control, parameter estimation, and otherareas. It has been extended to include nonlinear state models, nonlinearmeasurement models, non-white state and measurement noise, andcorrelation between state and measurement noise.

[0005] A succinct description of the Kalman filter can be found inApplied Optimal Estimation by A. Gelb (ed.) from The MIT Press, 1974 pp.1-179 which is herein incorporated by reference. A summary of the basicfilter and some notes on the notation follow from Gelb. The notation tobe used in the description is given in Table 1. A superscripted “T”after a matrix is the matrix transpose, and a superscripted “−1” after amatrix is it's inverse. Again following from Gelb, the discrete KFprocessing steps are given in Table 2. These steps advance the stateestimate and it's associated error covariance matrix from one time toanother. TABLE 1 Notation for Kalman Filter Description X_(t) ₁ _(|t) ₂State vector estimate for the system at time t₁ using measurements up toand including those from time t₂. Φ State transition matrix Q ProcessNoise Matrix R Measurement Noise Covariance Matrix P_(t) ₁ _(|t) ₂ Stateerror covariance matrix at time t₁ using measurements up to andincluding time t₂ Z Measurement vector, a set of vector measurements fora given time. H The measurement model matrix. This is the partialderivative of the measurement vector with respect to the state vector. IThe identity matrix

[0006] TABLE 2 Five Discrete Kalman Filtering Processing Steps Step StepName Step Processing 1 State Propagation X_(t+Δt|t) = ΦX_(t|t) 2Covariance Propagation P_(t+Δt|t) = ΦP_(t|t)Φ^(T) + Q 3 Gain CalculationK = P_(t+Δt|t)H^(T) (HP_(t+Δt|t)H^(T) + R)⁻¹ 4 State UpdateX_(t+Δt|t+Δt) = X_(t+Δt|t) + K(Z − HX_(t+Δt|t)) 5 Covariance UpdateP_(t+Δt|t+Δt) = (I − KH)P_(t+Δt|t)

[0007] An implementation issue with the KF is the matrix inverse calledfor in step 3, the gain calculation, from table 2. Computing the inverseis often a substantial burden, particularly in a real-time system Acommon practice, therefore, has been to divide up the processing of thevector measurement (steps 3, 4, and 5 from Table 2) into a sequence ofscalar measurement processing. If the measurement vector Z, has oneelement, in other words when we have a scalar measurement, themeasurement derivative with respect to the state, H, is a row vector. Inthis case the gain calculation, step 3 from Table 2 above reduces to adivision rather than a matrix inversion.$K = \frac{{PH}^{T}}{{HPH}^{T} + R}$

[0008] If a measurement vector, Z, contains multiple elements, it can bepartitioned into a set of scalar measurements. In the case of a diagonalmeasurement noise covariance matrix R, partitioning into scalarmeasurements is no more complicated than picking out the diagonal termfrom the matrix R, and processing it with the corresponding row of H andthe corresponding element of the vector Z. This is generally referred toas sequential measurement processing. We can expand our notation toinclude sequential measurement processing such that Z_(A), Z_(B), . . .are elements of the original vector measurement Z; H_(A), H_(B), . . .are the corresponding rows of H; R_(A), R_(B), . . . are thecorresponding diagonal elements of the noise covariance matrix. Thesequential measurement processing KF then follows these steps. The gainand residual calculations have been broken out for discussion of theChi-square test in subsequent paragraphs. The state and covariancenotation has been augmented for sequential processing, such thatX_(t+Δt|t,A) is the state vector estimate for time t+Δt givenmeasurements up to time t and measurement A. This state, updated formeasurement A, is the input state for steps 10 and 13 of measurement Bprocessing, and so on. TABLE 3 Discrete Sequential Kalman FilterProcessing Steps Repeated Step Steps Step Name Step Processing 1 StarePropagation X_(t+Δt|t) = Φ_(t|t) 2 Covanance Propagation P_(t+Δt|t) =ΦP_(t|t)Φ^(T) = Q Loop for each 3, 9, 15 . . . Residual Variance Γ_(A) =H_(A)PH_(A) ^(T) + R_(A) measurement: 4, 10, 16 . . . Residual δ_(A) =Z_(A) − H_(A)X_(t+Δt|t) A, B, C . . . 5, 11, 17 . . . Chi-Square teststatistic $\chi_{A}^{2} = \frac{\delta_{A}^{2}}{\Gamma_{A}}$

6, 12, 18 . . . Gain Calculation K_(A) = P_(t+Δt|t)H_(A) ^(T)/Γ_(A) 7,13, 19 . . . State Update X_(t+Δtt,A) = X_(t+Δt|t) + K_(A)δ_(A) 8, 14,20 . . . Covariance Update P_(t+Δt|t,A) = (I − KH)P_(t+Δt|t)

[0009] In the case where R is not diagonal, in other words themeasurement noises are correlated, techniques may be used whichtransform the R matrix to diagonal form for easy partitioning. Thenecessary transformations are defined by G. Bierman in FactorizationMethods for Discrete Sequential Estimation, Academic Press 1977, pp.47-55 which is herein incorporated by reference.

[0010] It is possible to prove that the order in which a set ofsequential measurements is processed in a linear system using the stepsof Table 3 is irrelevant to the final state estimate. A simpleone-dimensional example of this fact is given in Table 4. In Table 4,the two measurements are taken to be measurements of the state vectordirectly, thus H is a two-by-two identity matrix. The example is for thesequential measurement processing, not the state propagation, so itcorresponds to steps 3 through 14 inclusive of Table 3. The Chi-squaretest in the example of Table 4, will be described subsequently. Thereare two processing columns in Table 4, one for measurement order A,B,and the other for B,A. While the intermediate steps are different, thefinal state and covariance estimates are the same. Thus, in the priorart, the order in which measurements from the same time have beenprocessed in a sequential filter has been a matter of convenience, themeasurements are processed in whatever order they are gathered from themeasuring device or stored in memory, or some other mechanism. TABLE 4Single State KF Sequential Measurement Processing Example State Truth: X= 0 Estimate: X_(t|t−Δt) 32 1 State Variance P_(t|t−Δt) = 1 (stateestimate correctly modeled) Measurement Model H_(A) = H_(B) = H = 1R_(A) = R_(B) = R = 1 Measurements Z_(A) = −1 (One sigma measurement,correctly modeled) Z_(B) = 2 (Two sigma measurement, correctly modeled)Measurement Processing A, B B, A Order First Measurement Residual Γ_(A)= 1 + 1 = 2 Γ_(B) = 1 + 1 = 2 Variance First Measurement Residual δ_(A)= Z_(A) − X_(t|t−Δt) = −2 δ_(B) = Z_(B) − X_(t|t−Δt) = 1 FirstMeasurement Three (−2)² < (3² · 2) ? (1)² < (3² · 2) ? Sigma Chi-SquareTest Passed Passed First Measurement Gain K_(A) = P_(t|t−Δt)/Γ_(A) = 1/2K_(B) = P_(t|t−Δt)/Γ_(B) = 1/2 First Measurement State X_(t|t−Δt,A = X)_(t|t−Δt) + K_(A)δ_(A) = 0 X_(t|t−Δt,B = X) _(t|t−Δt) + K_(B)δ_(B) = 3/2Update First Measurement State P_(t|t−Δt,A) = (1 − K_(A))P_(t|t−Δt) =1/2 P_(t|t−Δt,B) = (1 − K_(B))P_(t|t−Δt) = 1/2 Variance Update SecondMeasurement Residual Variance$\Gamma_{B} = {{\frac{1}{2} + 1} = \frac{3}{2}}$

$\Gamma_{A} = {{\frac{1}{2} + 1} = \frac{3}{2}}$

Second Measurement δ_(B) = Z_(B) − X_(t|t−t,|A) = 2 δ_(A) = Z_(A) −X_(t|t−Δt,B) = −5/2 Residual Second Measurement Three Sigma Chi-SquareTest $\begin{matrix}{(2)^{2} < {\left( {3^{2} \cdot \frac{3}{2}} \right)?}} \\{Passed}\end{matrix}\quad$

$\begin{matrix}{\left( {- \frac{5}{2}} \right)^{2} < {\left( {3^{2} \cdot \frac{3}{2}} \right)?}} \\{Passed}\end{matrix}\quad$

Second Measurement Gain K_(B) = P_(t|t−Δt,A)/Γ_(B) = 1/3 K_(A) =P_(t|t−t,B)/Γ_(A) = 1/3 Second Measurement State X_(t|t) =X_(t|t−Δt,A) + K_(B)δ_(B) = 2/3 X_(t|t) = X_(t|t−Δt,B) + K_(A)δ_(A) =2/3 Update State Measurement State P_(t|t) = (1 − K_(B))P_(t|t−Δt,A) =1/3 P_(t|t) = (1 − K_(A))P_(t|t−Δt,A) = 1/3 Variance Update

[0011] One of the problems apparent to Kalman filter practitioners isthat measurements are often generated by physical mechanisms which canresult in measurement faults. These faults, by definition, arecomponents of the measurement values that are not included in themeasurement model. So, for example, if the measurement is a pseudorangemeasurement to a GPS satellite, and the direct ray to the satellite hasbeen blocked but a reflected ray is tracked, the resulting measurementhas a fault, errors that fall outside the model used to process themeasurement. This condition is shown in FIG. 2. Measurement faults haveat least two adverse affects; they corrupt the resultant estimate of thestate, and their effects are not considered in the update of the errorcovariance matrix. After processing a measurement with a fault, thestate error covariance matrix may no longer be an accuraterepresentation of the state error.

[0012] This problem led to the development of Chi-square measurementfault detection, sometimes referred to as Chi-square residual testing orinnovations variance testing. The purpose is to reject measurements thatwe believe have been corrupted by a fault, and therefore fall outside ofour filter model. The test is a statistical hypothesis test for a randomvariable. The application of a Chi-square test to the measurementresidual dates back, at least, to Mehra and Peschon [5] in a 1971 papertitled “An Innovations Approach to Fault Detection and Diagnosis inDynamic Systems” in Automatica, Vol. 7, pp. 637-640, 1971.

[0013] The chi-square residual test is a test on the hypothesis that themeasurement we are about to process is actually a member of themeasurement set we have defined. We define the residual as:

δ=Z−HX _(t+Δt|t)

[0014] The hypothesis is tested by forming the ratio of the observedresidual to it's expected value, thus, for a scalar measurement.$\chi^{2} = \frac{\delta^{2}}{{HPH}^{T} + R}$

[0015] The test would be passed if the computed value of chi-squared isless than a threshold value. A typical threshold of nine, for example,would be a three sigma test on the measurement, meaning we would expect99.7% of nor distributed measurements to pass the rest. So if the valueχ² computed above is greater than nine, we have detected a failure atthe three sigma level. This test is well-known and in common practice.

[0016] Since the residual variance, HPH^(T)+R, is also part of thescalar gain calculation, the chi-square test is often conducted in thesequential processing of the scalar measurements. Referring back toTable 3, the residual variance has been computed in step 3, and theChi-square test statistic was formulated in step 5. Thus the test iscomputed in-line with the sequential measurement processing. The examplein Table 4 shows the application of the test, the ratio is computed andchecked before the gain computation and the state update are completedfor each measurement. A flowchart of this process is given in FIG. 1.

[0017] The term Satellite Positioning System (SATPS) is used to define asatellite-based system transmitting ranging signals for use innavigation. Examples of SATPS include the Global Positioning System(GPS), the Russian GLONASS system, and the proposed European Galileosystem. In all of these systems the receiver develops pseudorangemeasurements to the transmitting satellites by synchronizing a codereplica with the transmitted satellite code. Doppler measurements, andintegrated doppler measurements, may be made as well, with integrateddoppler measurements commonly being referred to as delta-pseudorangemeasurements. Since the inception of the GPS, Kalman filtering has beena popular method for estimating the navigation state from thepseudorange and delta-pseudorange measurements. The paper “AircraftNavigation with the Limited Operational Phase of the NAVSTAR GlobalPositioning System” by L. R. Kruczynski, in Global Positioning System,Volume 1, pp. 154-165, 1980 from the Institute of Navigation, from his1977 conference paper, describes a Kalman filter used to process GPSmeasurements. In satellite positioning system measurement processingthere is some prior art on measurement fault detection methodology. Muchof this art is concerned with ad hoc attempts to define a measurementedit or weighting strategy to reduce the effects of multipath. In U.S.Pat. No. 5,883,595 by Colley, for example, a mechanism for de-weightingsuspect measurements is disclosed based on the size of their residualsand an undisclosed reliability factor. In U.S. Pat. No. 5,590,043 byMcBurney some adaptation of filter parameters on the basis of multipathpresence is disclosed, but with no mechanism apparent for detection ofthe multipath. The existence of these methods indicates the relativelack of success of unaltered Chi-square residual thong at least in thefield of satellite navigation, and calls out the need for improvedmeasurement detection and editing.

[0018] It can be seen then, that there is a need in the art for improvedmeasurement fault detection in state estimation. It can also be seenthat there is a need in the art for improved multipath detection insatellite positioning system receivers.

SUMMARY OF THE INVENTION

[0019] To minimize the limitations in the prior art, and to minimizeother limitations that will become apparent upon reading andunderstanding the present specification, the present invention disclosesa method and apparatus for improving measurement fault detection in asequential measurement processing estimator.

[0020] An object of the present invention is to provide for improvedmeasurement fault detection in state estimation. Another object of thepresent invention is to provide for improved multipath detection insatellite positioning receivers.

BRIEF DESCRIPTION OF THE DRAWINGS

[0021]FIG. 1 is a flowchart of a generic sequential measurementprocessing system for a Kalman filter with Chi-square residual testing;and

[0022]FIG. 2 is a diagram of a constellation of visible satellitepositioning system spacecraft being tracked by a terrestrial used with ameasurement corrupted by multipath in the so-called urban canyonenvironment.

DETAILED DESCRIPTION OF THE DRAWINGS

[0023] In the following description of the preferred embodiment,reference is made to the accompanying drawings which form a part hereof,and in which is shown by way of illusion a specific embodiment in whichthe invention may be practiced. It is to be understood that otherembodiments may be utilized and structural changes may be made withoutdeparting from the scope of the present invention.

[0024] Overview

[0025] The invention improves measurement fault detection by orderingthe processing of measurements in order of increasing probability of ameasurement fault, where we have some means to assess the relativelikelihood of a fault on the measurements to be processed, other thanthe measurement residual itself. For example, suppose a set ofpseudorange measurements is to be processed in a satellite positioningsystem receiver to obtain an estimate of the current position, as inFIG. 2. It appears in this Figure, and it is quite true that satellitesat a low elevation angle as seen from the receiver are much more likelyto have their transmitted signals affected by multipath and blockage fora terrestrial user. Thus we can say, without quantifying the exactprobability of failure, that the lower elevation satellite has a higherprobability of a multipath-induced fault on its pseudorange measurementthan a higher elevation satellite.

[0026] Similarly, multipath corrupted measurements are typicallyobserved to have a lower carrier to noise density ratio than uncorruptedmeasurements, and carrier tracking of a corrupted measurement may resultin larger phase-errors at the output of the carrier tracking loop. Thusthese two additional parameters, the C/N₀ of the measurement and thenumber and magnitude of the carrier tracking errors observed over somemeasurement interval are other indicators that may be used to ordermeasurements for sequential processing with increasing measurement faultprobability. These ordering elements could be used individually or incombination to determine the best processing order.

[0027] The key to the invention is to order the sequence of themeasurements presented to the filter for processing so that those with alower fault probability are processed first. This processing reduces thetrace of the state error covariance matrix and increases the probabilityof fault-detection on subsequent measurements. This has not beenrecognized in the prior art, Chi-square measurement fault detectionintroduced order dependence into sequential Kalman filter measurementprocessing. Once we have a measurement with a fault, the measurementprocessing sequence is no longer irrelevant. This is demonstrated inanother one-dimensional example in Table 5. Here we have a measurementwith a fault, and the outcome of the processing cycle is highlydependent on the order in which the measurements are processed. In thiscase, processing the measurements in the wrong order (B,A) results intwo estimation system failures, first, the faulty measurement is notdetected, then, following the processing of the faulty measurement, agood measurement is rejected because of the damage done to the stateestimate. TABLE 5 Single State KF Sequential Measurement Processing witha Measurement Fault Initial State Truth: X = 0 Estimate: X_(t|t−Δt) 32 1Initial State Variance P_(t|t−Δt) = 1 (state estimate error is correctlymodeled) Measurement Model H_(A) = H_(B) = H = 1 R_(A) = R_(B) = R = 1Measurements Z_(A) = −1 (One sigma measurement, correctly modeled) Z_(B)= 5 (Two sigma measurement − faulty) Measurement Processing A, B B, AOrder First Measurement Residual Γ_(A) = 1 + 1 = 2 Γ_(B) = 1 + 1 = 2Variance First Measurement Residual δ_(A) = Z_(A) − X_(t|t−Δt) = −2δ_(B) = Z_(B) − X_(t|t−Δt) = 4 First Measurement Three (−2)² < (3² · 2)? (4)² < (3² · 2) ? Sigma Chi-Square Test Passed Passed FirstMeasurement Gain K_(A) = P_(t|t−Δt)/Γ_(A) = 1/2 K_(B) = P_(t|t−Δt)/Γ_(B)= 1/2 First Measurement State X_(t|t−Δt,A = X) _(t|t−Δt) + K_(A)δ_(A) =0 X_(t|t−Δt,B = X) _(t|t−Δt) + K_(B)δ_(B) = 3 Update First MeasurementState P_(0|A) = (1 − K_(A))P₀ = 1/2 P_(0|B) = (1 − K_(B))P₀ = 1/2Variance Update Second Measurement Residual Variance$\Gamma_{B} = {{\frac{1}{2} + 1} = \frac{3}{2}}$

$\Gamma_{A} = {{\frac{1}{2} + 1} = \frac{3}{2}}$

Second Measurement δ_(B) = Z_(B) − X_(t|t−t,|A) = 5 δ_(A) = Z_(A) −X_(t|t−Δt,B) = −4 Residual Second Measurement Three Sigma Chi-SquareTest $\begin{matrix}{(5)^{2} < {\left( {3^{2} \cdot \frac{3}{2}} \right)?}} \\{Failed}\end{matrix}\quad$

$\begin{matrix}{\left( {- 4} \right)^{2} < {\left( {3^{2} \cdot \frac{3}{2}} \right)?}} \\{Failed}\end{matrix}\quad$

Updated State Estimate X_(t|t) = 0 X_(t|t) = 3 Second Estimate Error 0 3

[0028] An example of the invention has been given for sequentialmeasurement processing of measurements taken at a single time, however,the method can be extended. Suppose there are multiple measurementsavailable, at different times, but all available prior to the nextrequired reporting time for the state estimate. We can extend theconcept of ordered measurement processing to measurements taken atdifferent times by propagating the state and state error covarianceestimates backwards and forwards in time to process the measurementswith the lowest probability of a fault first. If, for example, we have ascalar measurement at time t₁ and a second measurement at a later timet₂, and we are not required to produce an updated state estimate until athird time, t₃. We can process these non-simultaneous measurements inorder by increasing fault probability even if the measurement at thelater time, t2, is the measurement with the lowest probability of afault.

[0029] The steps required to process the non-simultaneous measurementsare:

[0030] 1. Record measurement 1 at time t₁.

[0031] 2. Record measurement 2 at time t₂.

[0032] 3. Sort the measurements in order of increasing faultprobability, we'll assume for that the measurement with the lower faultprobability is the measurement 2.

[0033] 4. Propagate the state and covariance matrix forward from t₀ totime t₂

[0034] 5. Process measurement 2, updating the state and covariance attime t₂

[0035] 6. Propagate the state and covariance backwards in time to timet₁, noting that we subtract the relevant process noise matrix from thestate error covariance matrix when propagating backwards.

[0036] 7. Process measurement 1, updating the state and covariance attime t₁

[0037] 8. Propagate the state and covariance forward from t₁ to thereporting time t₃.

CONCLUSION

[0038] A method of improving measurement fault detection usingChi-square residual test by ordering measurements for sequentialprocessing by increasing probability of a measurement fault has beendisclosed. A simple one-dimensional example of the efficacy of themethod has been presented, and the application to ordering ofmeasurements for processing in a satellite positioning system receiverhas been presented.

[0039] In summary, the present invention provides a method for improvingfault detection. The method comprises recording a plurality ofmeasurements, each measurement having an associated measurement time,sorting the plurality of measurements, wherein the sort of the pluralityof measurements is performed in order of increasing fault probability,propagating a state matrix and a covariance matrix forward from thefirst time to a time of the measurement having a highest faultprobability, processing the measurement having the highest faultprobability first, continuing processing until all measurements havebeen processed, and propagating the state matrix and the covariancematrix forward from a last measurement time to a reporting time. Theprocessing comprises updating the state matrix and the covariance matrixat a time of the measurement having the highest fault probability, andpropagating the state matrix and the covariance matrix backwards in timeto a time of a measurement with a next highest fault probability.

[0040] The foregoing description of the preferred embodiment of theinvention has been presented for the purposes of illustration anddescription. It is not intended to be exhaustive or to limit theinvention to the precise form disclosed. Many modifications andvariations are possible in light of the above teaching. It is intendedthat the scope of the invention not be limited by this detaileddescription, but by the claims appended hereto.

What is claimed is:
 1. A method of fault detection, comprising:recording a plurality of measurements, each measurement having anassociated measurement time; sorting the plurality of measurements,wherein the sort of the plurality of measurements is performed in orderof increasing fault probability; propagating a state matrix and acovariance matrix forward from the first time to a time of themeasurement having a highest fault probability; processing themeasurement having the highest fault probability first, whereinprocessing comprises: updating the state matrix and the covariancematrix at a time of the measurement having the highest faultprobability; and propagating the state matrix and the covariance matrixbackwards in time to a time of a measurement with a next highest faultprobability; continuing processing until all measurements have beenprocessed; and propagating the state matrix and the covariance matrixforward from a last measurement time to a reporting time.
 2. The methodof claim 1 applied to a Global Positioning System (GPS) receiver.
 3. Themethod of claim 2, wherein the propagating the state matrix andcovariance matrix backwards in time further comprises subtracting aprocess noise matrix from a state error covariance matrix.
 4. The methodof claim 3, wherein the measurements are non-simultaneous.